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PREFACE 
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1.0  INTRODUCTION 

The  ambitious  operating  envelopes  desired  for  future  fighter  aircraft  dictates  an  aircraft 
industry  requirement  to  test  integrated  propulsion  system  components  at  a  ground  test  facility. 
In  response  to  this  industry  need,  the  Arnold  Engineering  Development  Center  (AEDQ  is 
developing  an  extensive  free-jet  test  capability  for  this  application  (Ref.  1).  One  factor  that 
inhibits  the  utilization  of  free-jet  testing  capability  at  the  AEDC  is  the  reliance  upon  empirical 
techniques  for  the  design  of  the  forebody  simulator  used  to  tailor  the  flow  field  to  the  desired 
flight  conditions.  To  alleviate  this  deficiency,  the  joint  ASD/AEDC  Aeropropulsion  System 
Test  Facility  (ASTF)  free-jet  development  technical  steering  committee  proposed  development 
of  an  aerodynamic  design  optimization  capability  applicable  to  the  forebody  simulator  design 
problem.  This  capability,  referred  to  as  the  "generation  6”  design  method  within  the  joint 
ASD/AEDC  ASTF  free-jet  development  plan,  was  intended  to  combine  nonlinear  optimization 
methods  with  the  powerful  analysis  capability  afforded  by  computational  fluid  dynamics 
(CFD).  Application  of  Euler  or  Navier-Stokes  CFD  analysis  codes  within  the  design 
optimization  is  motivated  by  the  aerodynamic  complexity  of  typical  ftee-jet  test  configurations. 
The  objective  of  this  report  is  to  document  the  progress  that  has  been  made  to  date  toward 
development  of  such  an  aerodynamic  design  optimization  capability. 

A  schematic  of  the  forebody  simulator  design  optimization  problem  is  illustrated  in  Fig. 
1.  This  figure  depicts  a  free-jet  inlet-engine  test  configuration  within  a  generic  ground  test 
facility  designed  to  evaluate  the  performance  of  an  integrated  propulsion  system.  The  design 
requirement  is  to  produce  a  flow  field  across  a  specified  inlet  reference  plane  in  the  free-jet 
installation  that  is  similar,  within  a  predefined  tolerance,  to  that  which  would  be  encountered 
in  flight.  One  proposed  way  of  achieving  this  is  by  appropriately  designing  a  ‘‘flow-tailoring" 
forebody  simulator,  and/or  varying  free-jet  flow  conditions  (total  pressure,  total  temperature, 
flow  angle,  and  Mach  number)  to  produce  the  desired  flow  field  at  the  inlet  reference  plane 
(Ref.  1).  The  plane  where  fluid  dynamic  similarity  is  required  will  hereafter  be  referred  to 
as  the  inlet  reference  plane  (1RP).  The  fundamental  assumption  of  this  free-jet  test  concept 
is  that  adequate  similitude  is  attained  whenever  the  fluid  dynamic  state  at  the  indicated  reference 
plane  adequately  matches  that  which  is  specified  for  a  given  operating  condition. 

The  task  of  designing  such  a  flow-tailoring  geometry  and  the  corresponding  test  conditions 
is  formidable.  The  test  designer  must  specify  a  forebody  simulator  geometry  and  free-jet 
fluid  properties  that  produce,  within  design  tolerance,  the  desired  fluid  dynamic  state  at  the 
region  of  interest.  Current  forebody  simulator  design  methods  rely  heavily  upon  prior 
experience  to  guide  a  trial-and-error  design  approach  using  subscale  testing  and  CFD  analyses 
to  evaluate  the  candidate  designs.  This  process  is  inefficient  and  does  not  ensure  an  orderly 
progression  toward  an  acceptable  design.  Thus,  the  specific  purpose  of  the  research  reported 
herein  was  to  develop  a  reliable  method,  based  upon  CFD  analyses,  for  the  specification 
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of  an  acceptable  set  of  aerodynamic  design  parameters  (both  geometric  and  fluid  dynamic) 
for  complex  designs  typical  of  those  encountered  in  an  aerodynamic  test  facility.  The  ultimate 
goal  of  the  research  is  to  optimize  the  design  of  a  complex  three-dimensional  (3-D) 
configuration,  such  as  illustrated  in  Fig.  I ,  in  a  timely  and  cost-effective  manner.  Development 
of  such  a  design  optimization  method  very  obviously  requires  the  ability  to  perform  an  accurate 
CFD  analysis  of  the  proposed  design.  However,  the  evaluation  of  various  CFD  techniques, 
relative  to  simulation  of  typical  free-jet  testing  configurations,  is  not  reported.  The  CFD 
evaluation  effort  is  being  conducted  in  parallel  with  the  design  optimization  development 
and  is  reported  in  AEDC-TR-90-21  by  M.  D.  McClure  and  J.  R.  Sirbaugh  (to  be  published). 

As  an  optimization  problem,  a  forebody  simulator  design  possesses  several  interesting 
features  in  that  (1)  the  flow-field  constraints  are  imposed  at  a  location  away  from  the  geometric 
surface  that  is  being  optimized;  (2)  both  fluid  dynamic  and  geometric  design  variables  must 
be  optimized;  (3)  discontinuous  or  localized  high-gradient  behavior  may  occur  within  the 
design  space  (e.g. ,  shocks  or  onset  of  separated  flow);  and  (4)  an  extensive  history  of  prior, 
similar  designs  does  not  exist.  Characteristics  (1)  and  (2)  do  not  constitute  a  well-posed  inverse 
design  problem;  however,  the  imposition  of  flow-field  data,  in  the  described  manner,  does 
properly  define  a  direct  optimization  problem.  Since  flow-field  data  are  prescribed  at  a  known 
reference  plane,  it  is  possible  to  define  an  objective  function  that  measures  a  norm  between 
the  reference  plane  flow  properties  associated  with  a  particular  design  point  (set  of  independent 
design  variable  values)  and  the  desired  reference  plane  flow  properties.  Since  the  norm  is 
a  function  of  the  given  design  variables,  the  optimization  task  is  to  determine  the  particular 
values  of  these  variables  that  produce  a  minimum  value  for  the  objective  function.  This  is 
accomplished  by  optimizing  the  selected  design  parameters  through  the  minimization  of  a 
nonlinear  least-squares  objective  function.  Flexibility  in  the  type  of  designs  that  can  be 
considered  is  provided  through  the  use  of  modern  CFD  codes  to  produce  the  function 
evaluations  required  by  the  optimization  algorithm.  The  resulting  direct  optimization  approach 
is  applicable  in  both  two  and  three  dimensions,  and  in  principle,  any  CFD  technique 
appropriate  to  the  flow  regime  of  interest  could  be  used. 

Using  an  Euler  or  Navier-Stokes  CFD  code  to  compute  design  space  gradients  within  an 
optimization  algorithm  has  received  little  prior  attention  in  the  literature.  It  is  demonstrated 
that  this  can  be  accomplished  by  applying  the  developed  design  technique  to  a  variety  of 
aerodynamic  design  problems.  The  test  problems  were  constructed  to  illustrate  the  applicability 
of  this  approach  to  realistic  designs  by  deliberately  selecting  poor  initial  conditions  for  the 
optimization  algorithm.  The  aerodynamic  optimization  examples  presented  include  a 
NACA0012  airfoil,  convergent/divergent  nozzles,  and  a  planar  supersonic  forebody  simulator. 
Although  the  method  applied  is  applicable  to  either  viscous  or  inviscid  flows,  only  one  viscous 
example  is  presented  because  of  the  increased  computational  expense  required  for  a  viscous 
CFD  analysis. 
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This  report  is  organized  so  that  the  design  optimization  technique  is  developed  and 
presented  in  Section  2.  The  design  technique  is  demonstrated  by  application  to  several 
aerodynamic  examples  in  Section  3.  Comments  relevant  to  the  application  of  this  technique, 
in  its  present  state  of  development,  are  presented  in  Section  4.  Lastly,  some  conclusions  relative 
to  design  optimization  using  the  developed  technique  and  some  recommendations  for  future 
research  are  presented  in  Sections  5  and  6. 

2.0  NUMERICAL  TECHNIQUE 

Using  CFD  to  optimize  aerodynamic  designs  is  currently  an  active  research  topic  in  the 
applied  mathematics  and  engineering  disciplines.  The  motivation  for  developing  and  using 
these  optimization  methods  in  the  design  process  is  to  reduce  the  overall  computational  effort 
needed  to  develop  aerodynamic  components  and  configurations,  which  will  optimize  a  selected 
measure  of  aerodynamic  performance.  Several  examples  of  aerodynamic  design  optimizations 
exist  in  recent  literature,  such  as  the  design  of  airfoils,  turbomachinery  cascades,  ducts,  and 
nozzles.  A  brief  literature  survey  was  provided  by  the  author  (Ref.  2)  and  is  not  repeated 
herein.  However,  it  is  noteworthy  that  similar  direct  optimization  methods  have  previously 
been  coupled  with  m eth od-of-characteristics  flow  solvers  at  the  AEDC  by  Varner  (Ref.  3) 
and  F.  L.  Shope  (unpublished  work).  Additionally,  a  similar  development  project  has  been 
proposed  recently  for  application  to  the  design  of  hypersonic  nozzle  contours  by  P.  F.  Hoffman 
(unpublished  work). 

2.1  BACKGROUND 

The  implementation  of  optimized  aerodynamic  design  generally  follows  one  of  three 
approaches:  (1)  inverse  design  methods,  (2)  basis  function  optimization  methods,  and  (3) 
direct  function  optimization  methods.  Generally,  true  inverse  design  methods  are  more  efficient 
than  either  the  basis  function  approach  or  direct  optimization  since  the  determination  of 
the  optimal  design  is  made  as  an  integral  part  of  the  CFD  analysis.  However,  since  many 
aerodynamic  design  problems  cannot  be  cast  in  an  inverse  form,  direct  optimization  methods 
and  basis  function  methods  are  often  applied. 

An  application  such  as  the  free-jet  forebody  simulator  illustrated  in  Fig.  1 ,  is  too  general 
for  successful  application  of  either  a  classic  inverse  method  or  the  basis  function  approach. 
However,  almost  any  design  problem  can  be  cast  as  a  direct  function  optimization  if  a  tangible 
measure  of  the  design's  quality  can  be  identified  to  define  an  objective  function  that  is 
responsive  to  changes  in  the  selected  design  parameters.  Thus,  the  implementation  of  a  direct 
aerodynamic  optimization  technique  is  investigated  and  reported  herein  as  well  as  in  Refs. 
2, 4,  and  5.  It  is  recognized  that  the  penalty  for  this  generality  is  a  potentially  less  efficient 
optimization  method  for  some  of  the  simpler  applications  that  may  be  of  interest  such  as 
airfoils  and  supersonic  nozzles. 
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With  current  computer  technology  and  CFD  algorithms,  many  complex  two-dimensional 
(2-D)  and  some  3-D  aerodynamic  designs  can  be  adequately  analyzed,  although  the 
computational  cost  can  be  very  high.  Applying  an  optimization  method  that  uses  CFD  to 
provide  function  evaluations  will  be  computationally  expensive  since  multiple  CFD  analyses 
are  required.  Even  so,  for  problems  such  as  the  previously  described  forebody  simulator, 
sonic  form  of  design  optimization  is  required  because  the  cost  of  the  available  alternatives 
(e.g.,  experimental  “cut  and  try”  in  a  wind  tunnel  using  an  inlet/forebody  simulator  model) 
may  be  even  more  prohibitive.  Additionally,  investigation  of  a  more  general  aerodynamic 
design  optimization  technique  will  help  prepare  the  way  for  future  enhancements  as  computer 
hardware,  CFD  algorithms,  and  optimization  algorithms  become  more  efficient. 

Within  this  research,  the  direct  optimization  problem  was  formulated  as  a  nonlinear  least- 
squares  minimization  using  existing  CFD  analysis  codes  to  provide  the  function  evaluations 
required  by  the  optimization  algorithm.  Both  Gauss-Newton  and  quasi-Newton  optimization 
algorithms  were  applied  to  minimize  the  least-squares  objective  function.  The  optimization 
algorithms  were  coupled  with  the  CFD  analysis  code  as  illustrated  in  Fig.  2  to  yield  the  desired 
interaction  between  the  CFD  analysis  capability  and  the  design  optimization  algorithm.  The 
optimization  code  was  kept  distinct  from  the  CFD  analysis  code  to  provide  the  analyst  with 
the  flexibility  to  select  the  most  appropriate  CFD  analysis  technique  for  a  given  design  problem. 

Implementation  of  this  optimization  technique  involved  three  primary  problems,  (1)  method 
of  function  evaluation,  (2)  selection  of  the  objective  function  and  implementation  of  the 
optimization  algorithm,  and  (3)  specification  of  the  design  parameters.  Each  of  these  items 
is  discussed  in  the  remainder  of  this  section. 

2.2  OBJECTIVE  FUNCTION  EVALUATION 

In  selecting  the  type  of  CFD  analysis  to  use  in  evaluating  the  objective  function,  the 
anticipated  flow  regime  to  be  encountered  computationally  was  identified  (Ref.  1).  A  typical 
free-jet  test  envelope  can  range  from  low  subsonic  flow  to  moderately  high  supersonic  flows, 
potentially  with  the  free-jet  nozzle  inclined  at  high  angles  of  attack  relative  to  the  test  article. 
The  appropriate  aerodynamic  analysis  for  the  motivating  problem  requires  a  complex,  3-D, 
flow-field  computation,  necessitating  the  application  of  an  Euler  code  or  a  Navier-Stokes 
code  to  produce  an  accurate  simulation.  However,  during  a  preliminary  design  phase,  less 
accurate  but  more  efficient  CFD  techniques  may  be  used.  By  formulating  the  optimization 
problem  as  a  nonlinear  least-squares  minimization,  the  particular  flow-field  analysis  technique 
applied  is  irrelevant  to  the  construction  of  the  optimization  algorithm  as  long  as  consistent 
and  repeatable  function  evaluations  are  obtained.  The  flow-field  simulation  must  be  consistent 
and  repeatable  in  the  sense  that  small  perturbations  to  design  parameters  are  accurately 
reflected  in  the  flow-field  solution.  This  is  important  because  the  implemented  optimization 
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algorithm  uses  these  function  evaluations  to  compute  design  space  gradients.  If  these  gradients 
are  inaccurate,  then  obviously  the  algorithm  would  not  converge  to  the  correct  solution. 

Although  the  direct  optimization  design  method  can  be  coupled  with  any  CFD  technique 
appropriate  to  the  problem  of  interest,  the  complexity  of  the  forebody  simulator  design 
problem  makes  it  necessary  to  demonstrate  that  the  direct  optimization  technique  can  be 
coupled  with  an  Euler  or  Navier-Stokes  solver.  Thus,  within  this  report  all  of  the  CFD  analyses 
were  made  using  PARC,  a  general  purpose,  finite  difference  Euler /Navier-Stokes  CFD  code 
(Ref.  6).  The  version  of  this  CFD  code  applicable  to  axisymmetric  and  2-D  configurations 
is  referred  to  as  PARC2D.  The  analogous  3-D  CFD  code  is  referred  to  as  PARC3D.  The 
PARC  codes  have  been  applied  at  the  AEDC  and  elsewhere  to  analyze  a  variety  of  complex 
internal  and  external  fluid  mechanics  problems  (Refs.  7  through  10).  This  particular  CFD 
code  was  selected  because  of  its  robustness,  ease  of  use,  and  reliability.  It  produces  consistent 
and  repeatable  flow  simulations  in  the  sense  that  small  perturbations  to  design  parameters 
are  accurately  reflected  in  the  flow  solution.  All  of  the  2-D  computational  grids  used  were 
generated  by  the  application  of  the  INGRID  code  developed  by  Soni  (Ref.  11). 

The  purpose  of  this  research  was  not  to  demonstrate  how  well  the  PARC  Euler/Navier- 
Stokes  code  can  simulate  a  particular  aerodynamic  phenomenon  or  to  improve  the  CFD 
analysis  capability  per  se,  but  to  demonstrate  that  an  Euler  /Navier-Stokes  code  can  be  coupled 
with  efficient  optimization  methods  to  produce  a  potentially  viable  technique  for  optimizing 
aerodynamic  designs  that  may  be  too  complex  to  design  by  other  available  means.  In  the 
aerodynamic  examples  presented,  no  effort  was  made  to  obtain  the  most  accurate  CFD 
simulation  for  the  given  configuration.  Whenever  possible,  a  minimum  number  of  grid  points 
were  applied  to  reduce  computation  time.  No  studies  were  made,  for  example,  to  assess  effects 
of  grid  distribution  on  the  CFD  simulation.  However,  it  was  necessary  to  monitor  the  level 
of  convergence,  particularly  at  the  reference  plane,  of  each  simulation  used.  It  was  necessary 
to  reduce  the  temporal  variation  of  flow  variables  at  the  reference  plane,  and  this  was 
accomplished  by  iterating  on  the  flow  field  until  the  norm  of  the  conservation  variables  at 
the  reference  plane  was  relatively  stationary,  typically  to  eight  significant  figures .  To  compute 
accurate  design  space  gradients  at  the  reference  plane,  the  dominant  change  in  the  residual 
must  be  attributable  to  the  change  in  the  design  parameters,  and  not  transient  effects. 
Demonstrating  that  CFD  can  be  used,  as  described,  within  an  optimization  algorithm  for 
aerodynamically  complex  configurations  will  extend  current  capabilities  in  aerodynamic  design 
optimization. 

2.3  OPTIMIZATION  ALGORITHM 

In  the  motivating  design  problem,  the  desired  fluid  dynamic  state  is  completely  known 
at  the  reference  plane,  either  from  experiment  or  free-stream  CFD  computation.  The  design 
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requirement  is  to  minimize  the  error  norm  between  these  quantities  and  corresponding  values 
computed  for  a  particular  design  point.  The  norm  was  selected  to  be  an  L4  norm  of  the 
difference  between  the  target  flow-field  variables  and  the  values  computed  for  a  given  design 
point.  The  nonlinear  least-squares  form  was  selected  because  (1)  the  method  is  very  versatile; 
(2)  extensive  literature  is  available  on  general  nonlinear  least-squares  minimization;  and  (3) 
efficient  Gauss-Newton  and  quasi-Newton  methods  are  well  documented  for  the  nonlinear 
least-squares  problem. 

The  nonlinear  least-squares  minimization  was  formulated  as  follows:  Let  the  residuals 
ri(P|,...,PM),  i=  1,2 . N,  be  functions  of  M  design  parameters, 

ri(Pi,P2 . Pm)  -  Yi  ~  fi(Pi,P2 . Pm)  (1) 

where  rj  denotes  the  difference  between  the  N  specified  reference  plane  quantities,  yj,  and 
the  corresponding  N  quantities  associated  with  the  M  parameters,  fj.  The  design  parameters 
may  be  geometric,  fluid  dynamic,  or  both.  To  minimize  rj,  in  the  least-squares  sense,  values 
for  the  parameters,  Pj,  are  found  that  minimize 

F(P„P2,...PM)  =  E  <ri(P|,P2,...,PM)}2  (2) 

i  =  1 


where  F  is  the  objective  function.  This  sum  can  be  written  in  vector  form  as  R(P)TR(P), 
where  P  denotes  a  vector  with  components  Pi  and  R(P)  denotes  a  vector  of  functions  with 
components  rj(P). 

Three  alternative  sets  of  fluid  dynamic  variables  were  considered  to  specify  the  reference 
plane  state  including,  (1)  the  Navier-Stokes  dependent  variables  in  conservation  form,  (2) 
the  Navier-Stokes  dependent  variables  in  nonconservation  form,  and  (3)  RP  total  conditions 
and  directional  Mach  number.  Preliminary  studies  conducted  in  this  research  using  a  simple 
airfoil  optimization  problem  detected  no  significant  difference  in  results  caused  by  the  choice 
of  dependent  variables. 


Except  when  otherwise  noted,  the  set  of  variables  used  herein  to  define  the  reference  plane 
(RP)  fluid  state  are  (1)  RP  total  pressure,  PTrp,  (2)  RP  total  temperature,  TTrp,  and  (3)  RP 
directional  Mach  number,  MIrp>  Myrp,  and  M^.  In  two  dimensions,  one  less  Mach  number 
component  is  required.  Variable  constraints,  when  applied,  are  imposed  by  adding  a  barrier 
function,  such  as  the  inverse  function  (Ref.  12)  to  the  objective  function.  Thus,  the  expression 
RTR  becomes 


Nr 

RTR  =  jEj  {(PTrp  -  Pt)?  +  (TTrp  -  Tt)2  +  (MXrp  -  MJ? 
+  (Myrp  "  M y)f  +  (M^  -  M*)?} 


(3) 
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where  Nr  is  the  number  of  reference  plane  points.  Since  each  term  in  Eg.  (3)  contains  five 
residual  fluid  dynamic  components,  in  order  to  put  this  in  the  form  of  Eq.  (2),  N  =  5Nr 
must  hold.  In  Eq.  (3)  the  subscript  rp  denotes  the  specified  reference  plane  values,  and 
unsubscripted  values  denote  reference  plane  values  computed  for  a  particular  trial  design 
(set  of  design  parameters).  Quantities  in  Eq.  (3)  are  normalized  by  appropriate  reference 
quantities  to  produce  target  reference  plane  values  of  order  one. 

A  popular  and  efficient  algorithm  for  minimizing  the  nonlinear  least-squares  form,  Eq. 
(1),  is  the  Gauss-Newton  method  (Ref.  12)  or  one  of  its  variants  such  as  Hartley’s  modified 
Gauss-Newton  method  (Ref.  13).  An  advantage  of  these  algorithms  as  applied  to  the  least- 
squares  form  is  the  elimination  of  the  need  for  the  Hessian  matrix  in  the  algorithm  formulation. 
Formation  of  the  Hessian  matrix  requires  specification  and  evaluation  of  N  x  M  x  (M 
+  l)/2  second  derivative  terms.  For  the  motivating  application,  computation  of  the  Hessian 
matrix  is  prohibitively  expensive  because  these  derivatives  must  be  approximated  by  finite 
differences.  Derivation  of  the  Gauss-Newton  method,  applied  to  the  least-squares  problem, 
is  available  from  several  sources  (Refs.  12  and  13)  but  is  repeated  in  Appendix  A. 

Applying  the  Gauss-Newton  method  to  minimize  Eq.  (1)  yields  an  optimization  algorithm 
of  the  form 


JTJAP  =  —  JTR 


(4) 


where  J  denotes  the  Jacobian  of  R  with  respect  to  P  defined  by 


(5) 


Equation  (4)  defines  an  M-by-M  system  of  equations  that  was  used  to  compute  the  change, 
AP,  in  the  design  parameter  solution  vector,  P.  To  apply  this  algorithm,  J  was  evaluated 
by  finite  difference  approximation  to  obtain  the  partial  derivative  of  each  residual  component 
with  respect  to  each  design  parameter.  This  requires  M  +  1  function  evaluations  to  compute 
the  M  partials  for  each  residual.  Since  a  CFD  solution  was  used  to  obtain  each  function 
evaluation,  approximation  of  this  Jacobian  was  by  far  the  most  expensive  part  of  the  algorithm. 


An  extension  of  this  algorithm  is  Broyden’s  quasi-Newton  method  (Refs.  12  and  14). 
Broyden’s  extension  modifies  the  standard  Gauss-Newton  method  by  approximating  the 
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Jacobian,  J,  at  the  k  +  1  iteration  strictly  from  the  Jacobian  and  other  data  available  at 
iteration  k  rather  than  recomputing  J  directly.  Since,  in  quasi-Newton  algorithms,  a  finite 
difference  approximation  to  J  is  made  only  for  the  initial  iteration,  the  accuracy  of  this 
approximation  is  even  more  important  than  for  the  Gauss-Newton  technique.  If  the  quasi- 
Newton  algorithm  is  to  converge  to  the  optimal  solution,  an  accurate  initial  approximation 
to  J  must  be  made.  Applying  Broyden’s  quasi-Newton  method  yields  an  optimization  algorithm 
identical  in  form  with  the  Gauss-Newton  method,  Eq.  (4),  and  is  given  by 


BT  BAP  =  -  BT  R 


(6) 


where  the  Jacobian  approximation,  Bt  at  iteration  k,  is  updated  for  iteration  k  +  1  according 
to 


Bk+1  =  Bk  + 


(ARk  -  BkAPk)  API 
(AP l  APk) 


(7) 


where  Bo  is  obtained  by  a  finite  difference  approximation  to  the  Jacobian.  This  modification 
to  the  Gauss-Newton  algorithm  is  based  upon  the  Jacobian  approximation  developed  by 
Broyden  (Ref.  14),  the  derivation  of  which  is  repeated  in  Appendix  B.  One  of  the  key 
assumptions  used  to  derive  Eq.  (7)  is  that  the  residual  change  in  directions  orthogonal  to 
the  direction  APk  predicted  by  Bk+i  is  identical  to  that  predicted  by  Bk.  The  imposition  of 
the  quasi-Newton  condition,  which  constrains  Bk+] to  hold  exact  derivative  information  in 
the  direction  of  APk  for  linear  R,  provides  the  other  conditions  necessary  to  uniquely 
determine  Bk  +  1.  Application  of  the  Gauss-Newton  algorithm  requires  M+l  function 
evaluations  for  each  iteration  since  the  Jacobian  is  approximated  by  finite  differences,  whereas 
Broyden 's  extension  requires  M  + 1  function  evaluations  for  the  first  iteration  but  only  one 
evaluation  for  subsequent  iterations.  Thus,  if  the  quasi-Newton  method  can  be  applied,  after 
the  first  iteration,  M  function  evaluations  are  eliminated  at  each  iteration.  For  typical 
aerodynamic  optimization  problems,  the  function  evaluation  is  the  dominant  part  of  the  cost, 
and  a  significant  savings  is  realized  whenever  the  quasi-Newton  algorithm  can  be  applied. 


It  is  well  known  that  the  Gauss-Newton  method,  using  analytical  derivatives,  determines 
a  search  direction  that  guarantees  a  reduction  in  the  objective  function  for  some  step  size 
in  that  search  direction  (Ref.  13).  Thus,  a  linear  search  technique  is  often  employed  once 
the  search  direction  is  determined  by  either  the  Gauss-Newton  or  quasi-Newton  algorithm. 
The  greatest  benefit  is  derived  from  the  linear  search  whenever  the  full  correction  computed 
by  the  optimization  algorithm  fails  to  produce  a  reduction  in  the  objective  function  value. 
Because  of  the  expense  of  function  evaluations,  a  comparatively  simple  method,  following 
Hartley  (Ref.  13),  was  applied  herein.  When  applying  Hartley’s  technique,  one  function 
evaluation,  in  addition  to  the  evaluation  at  the  predicted  optimum,  was  made  in  the  determined 
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search  direction.  The  base  design  point  and  the  two  new  design  points  along  the  search  direction 
were  then  used  to  define  a  quadratic  interpolation  function  that  was  analytically  solved  to 
yield  the  optimal  design  parameters  within  the  search  interval. 

The  linear  search  adds  at  least  one  additional  function  evaluation  per  iteration.  Because 
of  the  expense  of  function  evaluations,  the  linear  search  was  employed  only  when  a  full 
optimization  step  failed  to  produce  a  reduction  in  the  objective  function.  The  linear  search 
strategy  is  described  more  fully  in  Appendix  C. 

2.4  DESIGN  PARAMETERS 

For  the  problem  of  interest,  pertinent  design  variables  are  the  free-jet  fluid  dynamic 
parameters  and  the  variable  forebody  simulator  geometry.  The  free-jet  fluid  dynamic 
parameters  are  specified  as  jet  total  pressure,  jet  total  temperature,  and  jet  Mach  number. 
To  produce  an  efficient  optimization  method,  the  forebody  simulator  geometry  was  described 
parametrically  to  reduce  the  total  number  of  design  variables.  For  two-dimensional 
applications,  a  parametric  polynomial  representation  of  2-D  curves  as  given  by  the  Bemstein- 
Bezier  polynomial  (Ref.  15)  was  applied  as  follows: 

x(u)  =  E  — — —7-777-  ifll  ~  u)^-1)  xj;  0  s  u  s  1  (8) 

—  i  =  0  (n  -  i)!i! 

where  x  0,  x  i»  *2 . x_n  denote  the  position  vectors  of  the  n  + 1  geometric  control  points. 

In  this  form,  the  defined  curve  passes  identically  through  the  control  points  defined  by  vectors 
x  0  and  x  n  but  not  through  the  remaining  control  points.  This  allows  a  high  degree  of 
variability  for  a  given  number  of  design  parameters  relative  to  other  parametric  representations 
with  the  penalty  of  making  the  influence  of  each  parameter  upon  the  total  curve  somewhat 
obscure.  Application  of  the  Bezier  polynomials  prevents  the  large  oscillations  encountered 
when  using  interpolating  polynomials  because  of  the  "convex  hull”  property  of  Bezier 
polynomials.  This  property  ensures  that  the  Bezier  polynomial  lies  within  the  polygon  formed 
by  connecting  the  vertices  of  each  of  the  Bezier  control  vectors. 

For  3-D  applications,  an  extension  of  Eq.  (8)  defining  Bezier  surfaces  (Ref.  15)  was  applied. 
The  Bezier  surfaces  and  interior  were  defined  as  follows: 

x (u,v,w)  =  .  E  Q  j  Eft  E  0  gf  (u)  gj1  (v)  gk  (w)  Xjjk  (9) 

Here  x^k  denotes  the  position  vectors  of  the  control  points,  gf(u),  gjV),  and  gj,  (w)  are 
Bernstein  basis  functions  of  degree  p,q,  and  r,  respectively,  and  u,v,  and  w  are  parameters 
that  range  from  0  to  1.  The  Bernstein  basis  functions,  gf(u),  are  defined  by 
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(10) 


gf  (u)  =  — — Uj  (1  -  i  =  0,1,2,. ..p 

(P  —  1):I: 

with  the  other  basis  functions  analogously  defined. 

The  set  of  design  variables  are  of  diverse  type,  including  both  geometric  parameters  and 
fluid  dynamic  parameters.  Thus,  each  variable  was  nondimensionalized  by  an  appropriate 
reference  quantity.  The  nondimensionalization  was  chosen  to  make  each  design  variable 
nominally  of  order  one. 


3.0  NUMERICAL  EXAMPLES 

The  aerodynamic  design  optimization  technique  was  evaluated  by  optimizing  a  series  of 
numerical  examples  including,  (1)  algebraic  test  functions,  (2)  a  NACA0012  airfoil  in  inviscid 
flow,  (3)  a  NACA0012  airfoil  in  viscous  flow,  (4)  a  planar  inviscid  flow  in  a 
convergent/divergent  nozzle,  (S)  a  3-D  inviscid  flow  in  a  nozzle,  and  (6)  an  inviscid  supersonic 
flow  past  a  planar,  forebody  simulator.  In  each  of  these  examples,  the  optimization  problem 
was  formulated  so  that  a  global  minimum  exists  within  the  solution  space.  The  evaluation 
criterion  was  to  quantify  the  number  of  function  evaluations  required  to  determine  the 
optimum. 

The  primary  goal  of  this  research  was  to  demonstrate  the  feasibility  of  coupling  CFD 
analyses  capability  with  nonlinear  optimization  methods  to  produce  an  aerodynamic  design 
technique.  However,  if  the  developed  design  method  is  to  be  successfully  applied,  it  must 
also  be  efficient.  For  the  subject  application,  efficiency  can  be  equated  with  minimizing  the 
total  number  of  function  evaluations  required  to  isolate  the  optimum  since  this  is  where  the 
preponderance  of  the  computational  cost  is  incurred.  In  fact,  typical  computer  times  required 
to  evaluate  the  objective  function,  using  an  Euler  or  Navier-Stokes  solver  on  a  supercomputer, 
may  range  from  several  computer  minutes  for  a  simple  2-D  design  to  many  computer  hours 
for  a  complex  3-D  design  with  large  amounts  of  computer  memory  required.  The  optimization 
algorithm  applied,  which  used  the  function  evaluations,  typically  executed  in  5  to  10  sec  on 
a  CRAY*  XMP  supercomputer.  The  number  of  function  evaluations  was  minimized  by 
selecting  an  efficient  optimization  algorithm  and  applying  an  effective  geometric 
parameterization  to  reduce  the  number  of  design  variables. 

Madabhushi,  Levy,  and  Pincus  (Ref.  16)  developed  a  coupled  direct  optimization/CFD 
design  method  in  which  an  objective  function,  defined  as  the  average  duct  total  pressure 
loss,  was  minimized  as  a  general  nonlinear  function.  The  function  minimizations  were  made 
by  applying  the  Broyden,  Fletcher,  Goldfarb,  and  Shanno  (BFGS)  quasi-Newton  algorithm. 
The  selection  of  this  quasi-Newton  algorithm  was  based,  in  part,  upon  a  comparison  of  the 
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relative  performance  of  the  BFGS  algorithm,  a  conjugate  gradient  algorithm,  and  a  gradient 
algorithm  as  applied  to  the  minimization  of  ten  algebraic  test  functions.  The  relative 
performance  of  the  algorithms  was  measured  by  comparing  the  total  number  of  function 
evaluations  required  for  each  method  to  converge  to  the  global  minimum.  A  disadvantage 
of  the  BFGS  algorithm,  relative  to  forebody  simulator  optimization,  is  the  necessity  to  compute 
the  design  space  Hessian  matrix  (matrix  of  second  partial  derivatives).  Since  these  derivatives 
are  not  available  analytically,  computation  of  the  Hessian  represents  both  an  expensive 
computation  and  a  potentially  unreliable  computation  because  of  potential  inaccuracy  in  the 
objective  function  evaluation. 

Several  analytic  test  functions  were  used  by  the  author  (Ref.  2)  to  compare  the  relative 
efficiency  of  the  Gauss-Newton  algorithm  and  Broyden’s  quasi-Newton  algorithm,  as  applied 
to  the  nonlinear  least-squares  minimization  problem,  with  the  results  of  Madabhushi  et  al. 
The  goal  of  this  comparison  was  to  demonstrate  that  the  optimization  technique  applied  is 
reasonably  efficient  relative  to  other  available  optimization  algorithms.  If  the  comparisons 
are  favorable,  then  the  ease  of  application  and  versatility  afforded  by  the  nonlinear  least- 
squares  formulation  makes  this  an  attractive  technique  for  application  to  aerodynamic  design 
optimization.  Reiterating  a  further  advantage  of  the  least-squares  form  applied  herein  is  the 
elimination  of  the  computation  of  the  design  space  Hessian  matrix.  Generally,  the  results 
of  the  comparisons  were  favorable  with  the  advantage  of  utilization  of  a  less  complex 
optimization  algorithm.  Specific  details  are  given  in  Ref.  2. 

3.1  INVISCID  FLOW  ABOUT  A  PLANAR  AIRFOIL 

Optimization  of  a  planar  airfoil  in  inviscid  flow  by  specifying  "reference  plane"  (RP) 
values  at  a  station  downstream  of  a  NACA0012  airfoil  (Fig.  3)  provides  an  illustration  of 
the  design  optimization  technique.  The  RP  was  located  at  a  vertical  plane  beginning  at  the 
airfoil  trailing  edge  and  extending  five  chord  lengths  into  the  computational  domain.  The 
NACA0012  airfoil  is  a  well-known  airfoil  contour  that  has  been  extensively  analyzed  and 
is  defined  by 

y(x)  =  5t(0.2969x1/z  -  0.126x  -  0.35 16xz  +  0.2843x3  _  0.1015x4)  (11) 

where  the  parameter,  t,  specifies  the  airfoil  thickness.  For  the  NACA0012  airfoil,  the  thickness 
parameter  is  specified  as  0.12  for  a  chord  length  of  1.0089.  This  airfoil  contour  provided 
a  convenient  definition  of  a  one-parameter  design  optimization  problem  for  which  the  CFD 
analysis  was  very  simple.  The  purpose  of  this  example  was  to  provide  an  aerodynamically 
simple  design  problem  to  demonstrate  that  it  is  possible  to  optimize  an  aerodynamic  surface 
by  specifying  the  fluid  dynamic  state  at  an  RP  remote  to  that  surface. 
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The  PARC2D  CFD  code  was  used  to  define  the  target  RP  properties  by  computing  the 
inviscid  flow  Held  about  this  airfoil,  subject  to  the  boundary  conditions  indicated  in  Fig. 
3.  Nine  thousand  grid  points  were  used  to  resolve  the  domain.  Free-stream  properties  were 
held  constant  at  a  Mach  number  of  0.8  a  distance  of  five  to  ten  chord  lengths  away  from 
the  airfoil  surface.  Static  pressure  corresponding  to  the  specified  free-stream  Mach  number 
was  specified  at  the  indicated  computational  exit  plane.  The  RP  was  located  at  the  airfoil 
trailing  edge  and  extended  to  the  boundary  of  the  computational  domain.  The  influence  of 
the  body  was  shown  to  be  minima)  at  approximately  two  chord  lengths  into  the  domain. 
The  RP  properties  used  to  form  the  nonlinear  least-squares  objective  function,  as  defined 
by  Eq.  (2),  were  minimized  by  application  of  Broyden’s  algorithm. 

No  effort  was  made  to  obtain  a  highly  accurate  CFD  simulation.  However,  each  simulation 
was  scrutinized  to  assure  that  the  solution,  particularly  at  the  RP,  was  strongly  converged. 
In  this  example,  each  PARC2D  solution  was  converged  until  the  norm  of  the  RP  conservation 
variables  was  constant  to  eight  significant  figures.  For  this  simple  problem,  the  CFD  results 
were  consistent  with  those  obtained  by  Jameson  and  Mavriplis  (Ref,  17),  among  others,  as 
evidenced  by  the  pressure  coefficient  distribution  along  the  airfoil  surface  (Fig.  4). 

To  initialize  the  optimization,  the  airfoil  thickness  was  arbitrarily  perturbed  to  1/12  of 
its  original  value  as  an  initial  guess.  This  produced  a  very  flat  airfoil,  which  was  obviously 
distinct  from  the  target  profile  and  thus  produced  a  significant  perturbation  at  the  RP.  In 
fact,  as  seen  from  the  surface  pressure  profile  (Fig.  4),  the  target  airfoil  is  mildly  transonic 
with  a  weak  shock  appearing  at  approximately  midchord  on  the  airfoil  surface.  Conversely, 
the  almost  flat  profile  used  for  an  initial  guess  produced  very  little  distortion  of  the  free- 
stream  with  the  flow  field  remaining  subsonic  throughout  the  domain. 

The  optimal  design  parameter  was  derived  by  applying  Broyden’s  quasi-Newton  algorithm. 
For  the  first  iteration,  the  Jacobian,  J,  was  approximated  by  a  one-sided  finite  difference 
of  each  of  the  N  residuals.  This  difference  approximation  is  illustrated  in  the  following  example 
for  the  partial  derivative  of  the  i11*  residual  with  respect  to  the  design  parameter,  t, 

_3r|_  =  fj(t+At)  -  Tj(t) 

3t  At  1  ' 

Since  a  quasi-Newton  algorithm  was  applied,  a  good  approximation  to  the  design  space 
Jacobian  was  required  at  the  first  iteration.  This  Jacobian  was  not  recomputed  at  subsequent 
iterations  but  was  updated  approximately  according  to  Eq.  (7). 

To  give  an  indication  of  the  accuracy  of  this  Jacobian,  and  to  illustrate  that  stable 
computation  of  derivatives  was  possible  for  this  problem,  an  investigation  of  derivative 
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accuracy  versus  the  size  of  At  used  in  Eq.  (12)  was  performed.  Figure  5  shows  the  variation 
of  the  normalized  objective  function  derivative  with  respect  to  the  design  parameter  versus 
the  log  of  the  parameter  step  size  for  first-order  forward  differences  as  given  by  Eq.  (12). 
This  gradient  of  the  objective  function  was  used  to  provide  a  crude  assessment  of  the  partial 
derivatives  of  the  residual  components  that  were  used  in  the  optimization  algorithm.  Although 
this  does  not  provide  detailed  information  about  variation  of  individual  residual  derivatives, 
it  does  provide  a  means  of  measuring  global  variation  at  the  RP.  As  can  be  seen,  this  derivative 
was  sensitive  to  step  size,  but  approached  a  constant  value  for  step  sizes  less  than  0.001. 
In  fact,  the  actual  data  indicated  that  the  derivative  was  constant  within  0.S  of  1  percent 
for  step  sizes  less  than  0.001 .  Based  upon  this  analysis,  a  step  size  of  0.001  was  selected  for 
the  design  parameter. 

Figure  6  compares  the  target  geometric  profile  with  the  initial  guess  profile,  the  first  iteration 
profile,  and  the  optimal  profile  as  determined  by  Broyden’s  algorithm.  The  correct  geometry 
was  obtained  in  four  iterations,  which  required  five  function  evaluations  (CFD  solutions). 
Figure  7  shows  typical  variation  of  flow  variables  at  the  RP  as  evidenced  by  Mach  number 
profiles.  Figures  8  and  9  show  the  reduction  of  the  objective  function  and  the  convergence 
history  of  the  design  parameter,  t,  versus  iteration  number,  respectively.  As  evidenced  by 
these  figures,  the  Broyden’s  algorithm  isolated  the  global  minimum  quite  efficiently.  The 
data  from  which  Fig.  9  was  produced  indicates  that  the  optimum  was  located  within  1  percent 
in  2  iterations  and  was  isolated  within  0.1  of  1  percent  in  4  iterations. 

Airfoil  optimization  has  been  performed  by  several  other  researchers  by  prescribing  a 
pressure  distribution  along  the  airfoil  surface  and  solving  a  true  inverse  problem.  As  an 
interesting  example  of  the  versatility  of  the  nonlinear  least-squares  approach,  the  previously 
described  design  problem  was  also  solved  by  forming  the  least-squares  objective  function 
from  the  difference  between  the  airfoil  surface  pressures  for  a  given  design  point  and  a  specified 
pressure  distribution.  The  target  pressure  distribution  was  obtained  as  before  by  applying 
PARC2D  to  compute  the  inviscid  flow  field  about  this  airfoil,  subject  to  the  previously 
described  boundary  conditions  (Fig.  3).  The  airfoil  surface  pressure  distribution  was  then 
used  to  form  a  nonlinear  least-squares  objective  function,  defined  by 

RTR  =  (P,p  -  P)f  (13) 

where  P^  denotes  the  desired  airfoil  surface  pressure  at  one  of  the  N  RP  locations,  and  P 
denotes  the  airfoil  surface  pressure  for  a  given  design  point.  Optimization  of  the  airfoil  in 
this  manner  produced  results  very  similar  to  those  obtained  by  defining  a  trailing  edge  reference 
plane.  This  is  illustrated  in  Fig.  10  by  the  reduction  in  the  objective  function  and  in  Fig. 
11  by  the  convergence  of  the  design  parameter,  t.  For  this  formulation  the  optimum  was 
located  within  2  percent  in  4  iterations  and  was  isolated  within  0.3  of  1  percent  in  6  iterations. 
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3.2  VISCOUS  FLOW  ABOUT  A  PLANAR  AIRFOIL 

Formulating  the  aerodynamic  design  problem  as  a  nonlinear  least-squares  minimization 
allows  flexibility  in  selecting  the  type  of  CFD  simulation  as  well  as  allowing  application  to 
'  complex  designs.  For  example,  inclusion  of  viscous  effects  does  not  necessitate  any  changes 
to  the  applied  optimization  algorithm.  The  difference  in  the  overall  optimization  process 
occurs  only  in  the  CFD  analysis  step.  In  other  words,  a  more  elaborate  CFD  simulation 
provides  the  various  function  evaluations,  but  the  optimization  algorithm  does  not  recognize 
that  the  origin  of  the  function  evaluation  assumed  a  more  complex  physics  model. 

As  a  simple  demonstration  of  a  design  problem  in  which  viscous  effects  were  included, 
the  previously  described  NACA0012  airfoil  was  analyzed  with  function  evaluations  supplied 
by  viscous,  Navier-Stokes  simulations.  As  for  the  inviscid  flow  over  the  airfoil,  PARC2D 
was  used  to  define  the  target  RP  properties  by  computing  viscous  flow  about  a  NACA0012 
airfoil,  subject  to  the  indicated  boundary  conditions  (Fig.  12).  The  imposed  boundary 
conditions  were  identical  to  those  used  for  the  inviscid  airfoil  example,  except  the  airfoil 
surface  was  modeled  as  a  no-slip,  adiabatic  wall  rather  than  an  inviscid,  slip-wall  boundary. 
A  free-stream  Reynolds  number  of  106,  based  upon  chord  length,  was  specified,  which 
produced  an  attached  laminar  boundary  layer  (Fig.  13)  when  analyzed  with  PARC2D.  In 
this  example,  the  RP  was  extended  three  chord  lengths  into  the  domain  from  the  airfoil  trailing 
edge.  The  design  parameter  (airfoil  thickness)  was  perturbed  to  1/12  its  original  value  as 
an  initial  guess  to  begin  the  optimization.  Broy den’s  quasi-Newton  algorithm  was  then  applied 
to  derive  the  optimal  value  of  the  parameter.  Because  of  the  similarity  to  the  inviscid  problem, 
no  sensitivity  study  of  the  design  variable  was  performed. 

Figure  14  compares  the  target  geometric  profile  with  the  initial  guess  profile,  the  first 
iteration  profile,  and  the  optimal  profile  as  determined  by  Broyden’s  algorithm.  The  correct 
geometry  was  obtained  in  three  iterations,  which  required  five  function  evaluations  (CFD 
solutions).  Figures  15  and  16  show  the  reduction  of  the  objective  function  and  the  convergence 
history  of  the  design  parameter,  t,  to  the  known  optimum.  Again,  Broyden’s  algorithm  isolated 
the  global  minimum  very  efficiently.  The  data  from  which  Fig.  16  was  produced  indicates 
that  the  optimum  was  located  within  1  percent  in  2  iterations  and  was  isolated  within  0.1 
of  1  percent  in  3  iterations. 

This  viscous  example  was  modified  by  choosing  an  initial  guess  for  the  design  parameter 
that  was  double  the  correct  value,  producing  a  much  thicker  initial  airfoil  contour.  This  posed 
a  more  difficult  optimization  problem  because  the  computed  flow  Held  about  the  initial  guess 
airfoil  was  separated  (Fig.  17),  whereas  the  computed  flow  field  about  the  target  airfoil  was 
attached  for  the  prescribed  boundary  conditions  (Fig.  13).  The  absence  or  presence  of  separated 
flow  within  the  design  space  provided  an  abrupt  change  in  the  flow  field  for  candidate  designs 
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in  very  nearly  a  discontinuous  fashion.  To  assess  whether  the  optimization  technique  is  robust 
enough  to  optimize  designs  for  which  flow  separation  may  occur,  this  airfoil  was  analyzed 
with  the  RP  deliberately  placed  within  the  region  of  separated  flow.  The  RP  was  located 
at  the  trailing  edge  and  extended  from  the  airfoil  surface  to  the  edge  of  the  computational 
domain.  As  with  prior  examples,  PARC2D  was  used  to  provide  function  evaluations  at  the  RP. 

Convergence  for  this  example  was  similar  to  earlier  results,  although  somewhat  slower, 
with  the  optimum  design  variable  located  within  1  percent  in  5  iterations  and  to  within  0.1 
of  1  percent  in  6  iterations  requiring  7  function  evaluations.  Figure  18  illustrates  the  convergence 
to  the  target  geometry  by  comparing  the  target  geometric  profile  with  the  initial  guess  profile, 
the  first  iteration  profile,  and  the  optimal  profile  as  determined  by  Broyden’s  algorithm. 
Figures  19  and  20  illustrate  the  reduction  of  the  objective  function  and  the  convergence  of 
the  design  parameter,  t,  to  the  known  optimum,  respectively.  The  prominent  slope  change 
in  Fig.  19  coincides  with  the  first  candidate  design  for  which  minimal  separated  flow  was 
present. 

3.3  INVISCID  PLANAR  CONVERGENT/DIVERGENT  NOZZLE  FLOW 

Optimization  of  multiple  aerodynamic  design  parameters  was  illustrated  by  analyzing 
inviscid,  planar,  supersonic  flow  in  a  nozzle  (Fig.  21).  The  design  variables  for  this  problem 
were  inflow  total  pressure,  inflow  total  temperature,  and  the  nozzle  wall  contour  defined 
as  a  three-parameter  Bezier  curve.  The  vectors  used  to  define  the  Bezier  curve  were  located 
axially  at  the  inlet  plane,  the  midpoint,  and  the  exit  plane  as  indicated  in  Fig.  21.  The  exit 
plane  control  vector  was  held  constant  during  the  optimization.  The  y  coordinate  of  the  other 
two  vectors  were  allowed  to  vary  as  design  parameters  during  the  optimization  (Fig.  21). 
In  addition  to  these  two  geometric  parameters,  the  nozzle  inlet  total  temperature  and  the 
nozzle  inlet  total  pressure  were  allowed  to  vary,  yielding  a  total  of  four  design  parameters. 
Each  parameter  was  nondimensionalized  by  an  appropriate  reference  quantity.  The  Bezier 
parameters  were  normalized  by  the  target  nozzle  inlet  height,  total  pressure  by  the  target 
total  pressure,  and  total  temperature  by  the  target  total  temperature. 

The  simultaneous  variation  of  inlet  total  conditions  and  nozzle  wall  contour  is  not 
necessarily  representative  of  a  typical  nozzle  design  problem,  but  this  example  was  constructed 
because  simultaneous  variation  of  free-jet  total  conditions  and  a  variable  geometry  occurs 
with  the  supersonic  forebody  simulator  design  problem  that  motivated  the  present  research. 
The  target  solution  was  defined  by  selecting  the  Bezier  parameters  to  produce  a  nominal  2: 1 
area  ratio  nozzle  and  applying  PARC2D  to  solve  the  Euler  equations  using  1 ,600  grid  points 
to  resolve  the  domain.  To  form  an  initial  guess,  the  inflow  height  and  the  nozzle  throat  height 
were  reduced  so  that  the  initial  nozzle  area  ratio  was  approximately  doubled  and  the  nozzle 
throat  was  shifted  forward.  The  initial  guess  inflow  total  temperature  and  total  pressure  values 
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were  doubled.  The  disparity  between  the  target  design  and  the  initial  guess  is  illustrated  by 
comparing  the  centerline  Mach  number  profiles  for  the  target  design  and  for  the  initial  guess 
(Fig.  22).  As  can  be  seen,  the  exit  target  Mach  number  was  nominally  50  percent  below  the 
initial  guess  with  a  corresponding  variation  within  the  rest  of  the  nozzle.  Again,  no  special 
effort  was  made  to  produce  a  highly  accurate  simulation,  although  results  agree  well  with 
one-dimensional  theory  (Fig.  23). 


Broyden’s  quasi-Newton  algorithm  was  again  applied  in  this  example.  For  the  first  iteration, 
the  Jacobian  was  approximated  by  a  one-sided  finite  difference  of  each  of  the  N  residuals, 
as  illustrated  in  the  following  example  for  the  partial  derivative  of  the  1th  residual  with  respect 
to  the  j111  parameter 


A 

3Pi 


fj(Pl . Pj  +  APj,.-.,Pm)  fj(Pi,...,Pj,...,PM) 

APj 


(14) 


As  with  the  airfoil  example,  the  partial  derivative  of  the  objective  function  with  respect 
to  each  of  the  geometric  design  parameters  was  investigated  by  performing  a  sensitivity  analysis 
based  on  comparing  derivative  accuracy  versus  the  size  of  APj  used  in  Eq.  (14). 


Figure  24  depicts  the  variation  of  the  normalized  objective  function  differences,  as  given 
by  Eq.  (14),  with  respect  to  one  of  the  Bezier  design  parameters,  P|,  versus  the  log  of  the 
parameter  step  size.  As  can  be  seen,  this  derivative  is  sensitive  to  step  size,  but  attains  a  nearly 
constant  value  for  step  sizes  between  0.0001  and  0.002.  Between  these  limits,  the  partial 
derivative  is  constant  within  0.5  percent.  The  inaccuracy  observed  for  large  step  sizes  reflects 
the  nonlinearity  of  the  problem  and  the  inaccuracy  introduced  by  neglecting  higher  order 
terms  in  the  difference  approximation.  The  inaccuracy  for  very  small  step  sizes  results  from 
numerical  errors  in  the  function  evaluations  becoming  of  the  same  order  as  the  true  difference 
in  functional  value.  From  this  data  a  step  size  of  0.001  was  selected  for  the  Bezier  parameters. 
For  an  Euler  simulation,  the  residuals  vary  linearly  with  total  temperature  and  total  pressure; 
thus,  these  parameter  step  sizes  are  not  as  critical  in  this  example.  A  nondimensional  step 
size  of  0.001  was  selected  for  the  fluid  dynamic  parameters  to  be  consistent  with  the  geometric 
parameters. 


The  optimization  algorithm  converged  to  the  global  minimum  in  four  iterations,  which 
required  nine  function  evaluations.  Typical  convergence  of  the  reference  plane  properties 
is  illustrated  in  Fig.  25,  which  shows  reference  plane  Mach  number  profiles  for  the  target 
solution,  the  initial  guess,  the  first  iteration,  and  the  optimum.  The  reduction  in  the  objective 
function  is  illustrated  in  Fig.  26.  During  the  optimization,  inflow  total  conditions  converged 
to  within  1  percent  of  the  correct  value  in  1  iteration.  This  is  because,  for  this  example,  RP 
pressures  and  temperatures  scale  linearly  with  the  specified  total  conditions.  Thus,  the  residuals 
vary  linearly  with  these  parameters,  and  rapid  convergence  was  expected  since  a  Gauss-Newton- 
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type  algorithm  was  applied.  In  fact,  for  an  optimization  problem  in  which  the  residuals  are 
linear  in  all  of  the  parameters,  it  is  well  known  that  a  Gauss-Newton  algorithm  converges 
to  the  exact  solution  in  one  iteration.  The  geometric  parameters  converged  more  slowly  but 
still  at  an  acceptable  rate.  The  convergence  of  the  full  geometry  defined  by  the  individual 
parameters  is  illustrated  in  Fig.  27,  which  shows  the  iterative  variation  of  the  nozzle  wall 
contour.  The  convergence  of  the  individual  design  parameters  to  the  correct  solution  is 
illustrated  by  Figs.  28  through  31. 

3.4  IN  VISCID  THREE-DIMENSIONAL  NOZZLE  FLOW 

A  3-D  rectangular  nozzle  (Fig.  32)  was  used  to  demonstrate  that  the  nonlinear  least-squares 
optimization  method  is  applicable  in  three  dimensions.  The  nozzle  geometry  and  interior 
grid  were  defined  by  3-D  Bezier  polynomials,  Eq.  (9).  Four  control  points  were  specified 
at  each  of  five  axial  planes  so  that  each  axial  cross  section  was  rectangular.  The  entire  volume 
was  then  defined  using  the  coordinate  vectors  of  these  20  control  points  to  define  the  Bezier 
polynomials.  Two  design  parameters  were  defined  as  coefficients,  P|  and  P2,  that  controlled 
the  length  and  width  of  the  rectangular  cross  section  at  the  midplane  (Fig.  32).  The  target 
geometry  corresponded  to  values  of  unity  for  each  parameter,  which  produced  a  nozzle  with 
a  nominal  exit-to-throat  area  ratio  of  2.5.  Total  conditions  were  specified  at  the  nozzle  inlet. 
A  static  pressure  below  second  critical  was  selected  at  the  nozzle  exit.  This  allowed  fully 
supersonic  flow  to  develop  in  the  divergent  portion  of  the  nozzle.  These  geometry  and  boundary 
conditions  produced  a  flow  field  with  a  nominal  exit  Mach  number  of  2.5  when  analyzed 
with  the  Euler  version  of  the  PARC3D  code. 

For  an  initial  guess,  the  design  parameters  P}  and  P2  were  set  equal  to  2.0  and  2.5, 
respectively.  This  produced  a  nozzle  with  a  nominal  exit-to-throat  area  ratio  of  64.  A  PARC3D 
analysis  of  this  geometry  subject  to  the  described  boundary  conditions  produced  a  flow  field 
with  a  nominal  exit  Mach  number  of  5.8  using  23,000  grid  points  in  the  simulation.  Unlike 
the  target  nozzle,  which  was  square  at  each  axial  cross  section,  the  initial  guess  nozzle  had 
a  square  cross  section  at  the  inflow  plane,  which  transitioned  to  a  rectangular  cross  section 
at  the  midplane,  and  then  transitioned  again  to  a  square  at  the  exit  plane.  The  large  differences 
in  exit  flow  conditions  for  the  initial  guess  nozzle  compared  to  the  target  nozzle  were  selected 
to  illustrate  that  the  initial  guess  flow  Held  does  not  necessarily  need  to  closely  resemble  the 
desired  optimum  to  obtain  acceptable  results.  The  differences  in  the  flow  fields  for  the  target 
nozzle  and  the  initial  guess  nozzle  geometries  are  illustrated  by  comparing  the  centerline  Mach 
number  profiles  for  the  two  designs  (Fig.  33).  Also  presented  are  the  distributions  for  the 
first  iteration  and  for  the  computed  optimum. 

A  sensitivity  analysis  was  performed  on  design  parameter  Pi  that  indicated  that  parameter 
step  sizes  ranging  from  0.0001  to  0.01  produce  nearly  the  same  objective  function  derivative 
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(Fig.  34).  The  variation  in  Fig.  34  reflects  the  sensitivity  to  parameter  step  size  attributable 
to  nonlinear  effects  and  numerical  error  inherent  in  the  objective  function  evaluations.  A 
value  of  0.001  was  selected  for  the  design  variable  step  size  in  this  example.  Application  of 
Broyden’s  quasi-Newton  algorithm  did  not  converge  to  the  correct  optimum  for  this  test  case. 
It  is  suspected  that  the  Jacobian  computation  was  not  accurate  enough  for  reliable  application 
of  quasi-Newton  algorithm.  When  the  Gauss-Newton  algorithm  was  applied,  the  RP  properties 
converged  in  6  iterations  requiring  18  function  evaluations.  Figures  35  and  36  provide  an 
example  of  the  RP  convergence  by  comparing  y  and  z  centerline  Mach  number  profiles  at 
the  RP  for  the  target  solution,  the  initial  guess,  the  first  iteration,  and  the  computed  optimum. 
The  achieved  reduction  in  objective  function  and  the  design  parameter  convergence  is  depicted 
in  Figures  37  and  38,  respectively,  which  illustrate  that  the  global  minimum  of  the  objective 
function  was  isolated.  Convergence  is  further  demonstrated  by  comparing  the  wall  contours 
for  constant  y  planes  (Fig.  39)  and  constant  z  planes  (Fig.  40). 

3.5  INVISCID  SUPERSONIC  FLOW  ABOUT  A  PLANAR  FOREBODY  SIMULATOR 

A  3-D  analog  to  the  motivating  design  problem  was  constructed  and  is  shown  in  Fig. 
41.  The  configuration  shown  was  used  to  define  the  target  flow  variables  at  the  indicated 
reference  plane.  Thirteen  thousand  grid  points  were  used  to  resolve  the  domain,  which  was 
analyzed  inviscidly  by  the  application  of  PARC2D  in  the  Euler  mode.  In  this  example,  the 
Mach  number  at  the  inflow  plane  was  treated  as  a  design  variable  to  represent  the  variable 
free-jet  Mach  number  that  would  be  encountered  in  the  motivating  design  problem.  The 
forebody  simulator  geometry  was  defined  as  a  Bezier  curve  with  four  variable  parameters. 
The  geometric  design  variables  were  normalized  by  the  forebody  simulator  height,  and  the 
Mach  number  design  variable  was  normalized  by  the  target  Mach  number  to  make  each  design 
variable  of  the  same  order.  The  RP  was  located  one  inlet  height  in  front  of  the  inlet  entrance. 

The  target  design  was  a  relatively  blunt  forebody  with  an  incoming  Mach  number  of  3.0 
(Fig.  41).  The  reference  plane  was  located  downstream  of  the  detached  shock  emanating  from 
the  leading  edge  of  the  forebody.  Also  within  the  flow  field,  weak  shocks  reflect  from  the 
cell  wall  boundary  and  from  the  nozzle  walls.  The  simulation  was  run  with  an  exit  pressure 
low  enough  to  maintain  fully  supersonic  flow  across  the  computational  exit  plane.  The  five 
design  parameters  shown  in  Fig.  41  were  then  perturbed  to  initialize  the  optimization.  For 
the  initial  guesses,  the  Mach  number  was  increased  by  50  percent,  and  the  geometric  parameters 
were  reduced  by  10  to  70  percent  to  intentionally  provide  a  poor  initial  guess. 

The  disparity  between  the  target  design  and  the  initial  guess  configuration  is  illustrated 
by  comparison  of  the  respective  Mach  number  contours  (Fig.  42  and  43).  As  can  be  seen, 
the  shock  structures  in  front  of  the  reference  plane  are  very  different,  leading  to  large 
discrepancies  in  RP  flow  variables.  The  forebody  shock  passes  very  close  to  the  RP  in  the 
initial  guess  configuration,  which  further  complicates  the  optimization  task. 
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The  RP  objective  function  was  minimized  by  the  application  of  the  Gauss-Newton 
algorithm.  Analysis  of  the  objective  function  gradients  indicated  that  the  minimum  derivative 
variation  was  about  2  percent  (Fig.  44)  for  nondimensional,  geometric  variable  step  sizes 
between  0.02  and  0.05.  These  derivatives  were  not  accurate  enough  for  successful  application 
of  the  quasi-Newton  algorithm.  However,  when  the  Gauss-Newton  algorithm  was  applied, 
convergence  was  obtained  in  five  iterations  requiring  31  function  evaluations.  The  maximum 
Mach  number  deviation,  at  the  RP,  was  less  than  1  percent  after  only  two  optimization  steps 
requiring  13  function  evaluations,  which  for  many  applications,  may  be  adequate.  A 
comparison  of  the  reference  plane  Mach  number  profiles  for  the  initial  guess,  the  target 
solution,  and  the  final  converged  solution  is  made  in  Fig.  45,  which  shows  excellent  agreement 
between  the  target  and  final  solutions.  Figure  46  plots  the  objective  function  versus  design 
iteration  number  illustrating  that  the  global  minimum  of  the  objective  function  was  isolated. 
The  convergence  history  for  each  of  the  five  design  parameters  is  illustrated  in  Fig.  47. 

3.6  INVISCID  SUPERSONIC  FLOW  ABOUT  A  REDUCED  LENGTH  FOREBODY 

SIMULATOR 

The  motivating  design  problem  requires  that  RP  properties  be  matched  within  design 
tolerance  while  significantly  shortening  the  actual  aircraft  forebody.  An  example  emulating 
this  type  of  optimization  problem  was  constructed  by  seeking  a  forebody  half  the  length  of 
the  described  target  forebody  as  depicted  in  Fig.  41  and  retain  the  RP  target  values  defined 
in  Section  3.5.  This  length  restriction  excluded  an  exact  match  of  RP  variables  from  the 
optimization  design  space.  Each  candidate  design  was  analyzed  as  a  full  test  cell  configuration 
without  using  a  plane  of  symmetry  as  with  the  target  configuration.  Each  design  was  analyzed 
inviscidly  by  the  application  of  the  blocked  version  of  PARC2D  in  the  Euler  mode  using 
26,000  grid  points  to  resolve  the  domain  in  each  CFD  simulation. 

The  optimization  problem  was  constructed  using  four  design  variables.  The  inflow  plane 
Mach  number  (normalized  by  M  =  3.0  as  a  reference  Mach  number)  was  used  as  a  fluid 
dynamic  design  variable.  The  forebody  simulator  geometry  was  again  defined  as  a  Bezier 
curve  with  the  vertical  component  of  three  of  the  Bezier  control  vectors  used  as  design  variables. 
These  geometric  design  variables  were  normalized  by  a  reference  length  equal  to  the  height 
of  the  target  forebody.  As  with  the  full-length  forebody  design,  the  Gauss-Newton  algorithm 
was  applied.  In  three  iterations,  requiring  16  function  evaluations,  the  deviation  from 
normalized  RP  objective  function  components  of  total  pressure  and  axial  Mach  number  was 
as  indicated  in  Figs.  48  and  49.  In  these  figures,  each  RP  component  was  normalized  by 
the  corresponding  target  value  to  facilitate  comparison  of  the  various  design  iterations.  The 
target  values  for  the  vertical  Mach  number  component  were  very  near  zero;  thus,  comparison 
of  the  actual  residual  components  is  most  meaningful  and  is  presented  in  Fig.  50.  The  total 
temperature  component  of  the  objective  function  is  not  presented  since  total  temperature 
remains,  constant,  within  numerical  accuracy,  for  this  inviscid  computation. 
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As  seen  from  Figs.  48  through  SO,  each  component  of  the  objective  function  was  reduced 
from  its  initial  value.  The  convergence  and  rate  of  reduction  of  the  total  objective  function 
is  further  illustrated  in  Fig.  SI .  Further  iterations  failed  to  yield  any  appreciable  reduction 
in  the  objective  function.  Additionally,  the  contribution  from  the  total  temperature  components 
to  the  objective  function  sum  became  of  the  same  order  as  the  total  pressure  and  Mach  number 
components.  This  indicated  that  further  reduction  in  the  objective  function  was  unlikely, 
since  for  the  Euler  simulations,  any  mismatch  in  total  temperature  could  be  attributed  to 
numerical  error.  The  comparison  of  the  RP  Mach  number  for  various  iterations,  as  shown 
in  prior  examples,  is  presented  in  Fig.  52.  The  variation  of  the  geometric  contour  is  illustrated 
in  Fig.  S3.  In  this  two-dimensional  example,  the  Gauss-Newton  algorithm  was  able  to 
substantially  reduce  the  deviation  of  RP  fluid  dynamic  parameters  from  the  prescribed 
distributions  leading  to  a  substantially  improved  design.  The  number  of  required  function 
evaluations  is  also  within  practical  limits  using  current  technology. 

4.0  USAGE  CONCEPTS 

The  computer  program  that  was  used  to  generate  the  previous  results  is  available  for 
application  at  the  AEDC.  Although  it  is  still  a  developmental  code,  it  is  relatively  easy  to 
apply.  This  section  is  not  intended  to  serve  as  a  complete  user’s  guide  for  the  developmental 
code,  but  will  serve  to  briefly  describe  information  required  to  apply  the  design  optimization 
code  in  its  present  form. 

Input  variables  required  by  the  code  are  supplied  through  a  namelist  format.  However, 
before  defining  the  necessary  input  variables,  an  overall  description  of  the  design  process 
is  helpful.  Examination  of  Fig.  2  reveals  four  key  steps  in  the  design  optimization  process, 
(I)  definition  of  design  variables,  (2)  analysis  of  the  candidate  design,  (3)  measure  of  the 
design  quality,  and  (4)  adjustment  of  the  design  variables. 

Design  optimization  methods  are  intended  to  assist  the  engineer  in  developing  a  quality 
design.  They  are  not  intended  to  eliminate  the  need  for  good  engineering  judgment.  Defining 
the  design  variables  and  measuring  the  design  quality  are  two  areas  where  a  thorough 
understanding  of  the  design  requirements  and  a  recognition  of  the  important  design  parameters 
is  essential  to  developing  a  quality  design.  Within  this  research,  a  least-squares  norm  at  the 
IRP  was  used  as  an  objective  function  to  measure  the  design  quality.  This  objective  function 
was  assumed  to  be  a  function  of  the  selected  design  variables.  If  the  IRP  objective  function 
is  not  responsive  to  variation  of  the  selected  design  variables,  then  the  problem  is  not  well 
defined.  It  must  also  be  remembered  that  design  variable  specification  defines  the  allowable 
design  space  that  will  be  searched  for  the  optimal  design.  If  the  defined  design  space  is  too 
restrictive,  the  resulting  design  may  be  unsatisfactory  even  though  the  optimal  value,  within 
the  allowable  design  space,  was  properly  located.  This  is  one  of  the  reasons  that  the  Bezier 
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surfaces  were  selected  for  geometric  parameterization.  They  exhibit  the  desirable  features 
of  providing  a  wide  range  of  geometric  flexibility  for  a  relatively  small  number  of  design 
variables.  An  increase  in  the  number  of  design  variables  directly  increases  the  number  of 
function  evaluations  (CFD  simulations)  required  to  attain  an  optimum. 

As  with  any  CFD  simulation,  all  results  should  be  carefully  analyzed.  For  the  coupled 
optimization  method,  multiple  CFD  simulations  are  combined  to  determine  design  space 
gradients  that  iteratively  leads  to  an  optimized  design.  If  the  results  of  the  simulations  are 
not  consistent  and  reliable,  then  poor  results  are  to  be  expected  (garbage-in  yields  garbage- 
out).  Some  items  to  pay  particular  attention  to  include  (1)  large  variations  in  grid  quality 
between  simulations,  (2)  temporal  variations  attributable  to  inherent  unsteadiness  or  poor 
convergence,  and  (3)  generally  poor  results  near  the  IRP. 

For  an  N-dimensional  design  problem,  a  CFD  simulation  is  required  for  the  base  design 
and  N  perturbations  to  the  base  design.  The  perturbation  value  is  specified  by  the  user  and 
should  be  consistent  with  the  design  variable.  Typically,  this  value  is  determined  by  examining 
the  design  space  as  described  in  Section  3  of  this  report.  The  optimization  code  is  executed 
as  a  postprocessor  after  each  CFD  simulation  to  evaluate  the  objective  function  and  to  compute 
the  value  of  each  design  variable  to  be  used  for  the  next  function  evaluation.  One  initialization 
run  is  required  to  establish  target  IRP  values.  In  this  manner,  the  optimization  code  guides 
the  engineer  to  the  optimal  design. 

4.1  FILE  USAGE 

The  job  control  file  must  fetch  necessary  input  files,  compile,  load,  and  execute  the  source 
code  and  dispose  generated  output  files.  The  source  code  requires  four  input  files  as  follows: 

1 .  Unit  #2  —  a  PARC  format  restart  file  corresponding  to  the  current  set  of  design 
variables. 

2.  Unit  #3  —  the  optimization  code  restart  file,  which  contains  IRP  flow  variables 
and  IRP  grid  coordinates  for  all  prior  design  iterations. 

3.  Unit  #  5  —  a  namelist  input  file  that  contains  code  variable  definitions  to  control 
program  execution. 

4.  Unit  #9  —  a  design  variable  history  file  that  is  used  to  monitor  the  iterative 
behavior  of  the  design  variables  and  convergence  toward  the  optimal  design. 
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Two  output  files  are  created  containing 

1.  Unit  #7  —  the  updated  optimization  code  restart  file. 

2.  Unit  #  10  —  the  updated  optimization  code  design  variable  history  file. 

Coupling  the  optimization  code  to  an  analysis  code  other  than  PARC  would  necessitate  minor 
changes  in  the  input  subroutine. 

4.2  NONDIMENSIONALIZATION 

Since  diverse  flow  properties  are  used  to  comprise  the  IRP  residuals,  each  residual 
component  is  normalized  by  a  user-selected  reference  value,  as  noted  in  Section  2.3.2.  A 
general  recommendation  is  to  make  each  residual  term  of  order  one,  in  which  case  each 
component  of  the  residual  will  be  driven  to  zero  at  approximately  the  same  rate.  However, 
there  may  be  instances  where  the  user  wishes  to  bias  the  solution  by  weighting  one  set  of 
residuals  more  heavily  than  other  residuals.  This  can  be  accomplished  by  choosing  an 
appropriate  set  of  reference  properties. 

4.3  NAMELIST  CONTROL 

All  user  inputs  to  the  design  code  are  provided  from  Unit  #  5  through  the  namelist  titled 
“Control.”  These  inputs  are  used  to  define  (1)  PARC  nondimensionalization  variables,  (2) 
design  code  nondimensionalization  variables,  (3)  IRP  target  value  definition,  (4)  variable 
geometry  indices,  and  (S)  optimization  code  operation.  Some  variables  contained  within  this 
namelist  are  strictly  for  code  development  and  diagnostic  use.  Description  of  variables  needed 
to  execute  the  design  optimization  code  are 


DTH 

Real  vector  specifying  the  perturbation  value  for  each  design  variable 
during  derivative  computations 

1BACKUP 

Development  and  diagnostic  use 

IC1 

Integer  vector  identifying  the  set  of  candidate  designs  to  use  within 
the  current  optimization  step 

FMACH 

The  reference  Mach  number  used  to  normalize  the  directional  Mach 
number  residual  components 

GAMMA 

PARC  code  specific  heat  ratio 
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IBPATH 

i 

JEGEO 

JIRP 

JSGEO 

KEGEO 

KSGEO 

KEIRP 

KSIRP 

KSKIP 

LEGEO 

LSGEO 

LEIRP 

LSIRP 

NCASES 

NCYCBACK 


Integer  vector  specifying  the  set  of  candidate  designs  needed  to  form 
quasi-Newton  Jacobian  algorithm  (Broyden’s  update  formula) 

Beginning  j-index  for  the  variable  design  geometry.  Used  only  for 
postprocessing  of  results. 

IRP  j-index 

Variable  design  geometry  ending  j-index.  Used  only  for  postprocessing 
of  results. 

Variable  design  geometry  beginning  k-index.  Used  only  for 

postprocessing  of  results. 

Variable  design  geometry  ending  k-index.  Used  only  for 

postprocessing  of  results. 

IRP  ending  k-index 

IRP  starting  k-index 

Skip  factor  to  specify  the  frequency  with  which  a  full  Gauss-Newton 
optimization  step  is  performed.  KSKIP  =  1  implies  all  steps  will  be 
Gauss-Newton,  KSKIP  =  5  implies  every  fifth  step  is  Gauss-Newton 
with  all  other  steps  quasi-Newton,  etc. 

Variable  design  geometry  beginning  1-index.  Used  only  for 

postprocessing  of  results. 

Variable  design  geometry  ending  1-index.  Used  only  for  postprocessing 
of  results. 

IRP  ending  1-index 

IRP  starting  1-index 

The  cumulative  number  of  candidate  designs 
Development  and  diagnostic  use 
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NBROYSTP 

The  number  of  designs  used  to  execute  the  Broyden  update  formula 

NSTART 

Initialization  parameter.  If  NSTART =0  the  code  performs  a  first 
pass  initialization  to  define  IRP  values 

PREF 

PARC  code  reference  pressure 

PT  NORMAL 

The  reference  pressure  used  to  normalize  the  total  pressure  residual 
components 

TREF 

PARC  code  reference  temperature 

TT  NORMAL 

The  reference  temperature  used  to  normalize  the  total  temperature 
residual  components 

5.0  SUMMARY  AND  CONCLUSIONS 

A  direct  aerodynamic  design  optimization  technique  that  couples  an  existing  Euler/Navier- 
Stokes  solver  with  efficient  nonlinear  least-squares  minimization  algorithms  has  been  developed 
and  demonstrated  by  successfully  applying  the  technique  to  representative  aerodynamic  design 
problems.  This  demonstrated  that  design  space  gradients,  required  by  the  Gauss-Newton- 
based  optimization  algorithms,  could  be  reliably  formed  from  Euler/Navier-Stokes  CFD 
simulations.  In  fact,  for  two  examples  presented,  a  quasi-Newton  algorithm  was  successfully 
applied  to  determine  the  optimal  design  variables.  The  examples  were  deliberately  started 
from  poor  initial  designs  to  assess  the  capability  of  the  design  technique  to  converge  from 
poor  initial  conditions.  Although  this  does  not  yet  achieve  the  desired  3-D  forebody  simulator 
design  capability,  it  does  represent  significant  progress  toward  satisfaction  of  that  goal. 

To  date,  this  research  has  evaluated  the  feasibility  of  coupling  nonlinear  optimization 
methods  with  CFD.  Existing  Gauss-Newton  and  quasi-Newton  optimization  algorithms  were 
employed  with  minimal  modifications  with  good  results  obtained  for  several  representative 
aerodynamic  design  problems.  The  quasi-Newton  algorithm  was  not  successful  for  the  more 
aerodynamically  complex  examples,  but  was  significantly  more  efficient  when  applicable. 
This  suggests  that  alternating  between  the  two  algorithms  may  produce  a  more  efficient 
optimization  strategy  if  reliable  switching  criteria  can  be  determined.  Because  of  the  complexity 
of  the  motivating  application,  the  optimization  algorithms  were  coupled  with  an  Euler/Navier- 
Stokes  CFD  code,  which  makes  the  function  evaluations  computationally  expensive.  However, 
the  same  method  could  be  coupled  with  a  less  complex  CFD  technique  for  design  problems 
in  which  an  Euler  or  Navier-Stokes  simulation  is  not  required. 
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The  approach  presented  provides  the  designer  with  a  potentially  powerful  tool  to  assist 
in  many  designs  for  which  a  measure  of  design  quality  (objective  function)  can  be  adequately 
defined.  Because  of  the  flexibility  afforded  by  the  current  CFD  codes,  this  technique  can 
be  applied  to  virtually  any  steady-state  aerodynamic  optimization  problem  for  which  the 
selected  CFD  code  is  capable  of  providing  a  reliable  aerodynamic  analysis.  The  design  method 
is  reliable  in  the  sense  that  improvements  to  the  initial  design  are  achieved,  and  it  is  efficient 
in  terms  of  minimizing  the  number  of  function  evaluations  required  to  determine  the  optimal 
design  variables.  For  complex  2-D  aerodynamic  designs,  the  method  requires  on  the  order 
of  10  to  50  function  evaluations,  depending  upon  the  quality  of  the  initial  design,  number 
of  optimization  parameters,  and  RP  simulation  tolerances.  This  number  of  function 
evaluations  is  feasible  with  access  to  modern  supercomputers  and  use  of  appropriate  CFD 
codes  that  efficiently  yield  the  desired  flow  physics.  As  demonstrated  by  one  simple  example, 
the  design  method  is  applicable  in  three  dimensions.  However,  the  efficiency  of  the  optimization 
process  must  be  enhanced  before  this  technology  can  be  considered  a  practical  3-D  design  tool. 

6.0  RECOMMENDATIONS 

Achieving  the  “generation  6”  forebody  simulator  design  capability,  proposed  by  the 
ASD/AEDC  ASTF  free-jet  development  technical  steering  committee,  requires  continued 
development  of  aerodynamic  design  optimization  techniques.  Computational  fluid  dynamics 
is  certainly  not  a  mature  science  and  the  use  of  CFD  within  a  design  optimization  framework 
is  in  its  infancy.  If  continued  development  of  a  “generation  6“  design  capability  is  deemed 
a  priority  item  at  the  AEDC,  then  the  following  topics  need  to  be  investigated. 

1 .  Investigate  the  application  of  design  optimization  via  control  theory  as  proposed 
by  Jameson  (Ref.  18). 

2.  Investigate  the  application  of  genetic  algorithms  (Ref.  19)  to  minimize  the  number 
of  required  CFD  simulations. 

3.  Improve  the  present  direct  optimization  method  by  combining  the  efficiency 
characteristics  of  the  quasi-Newton  algorithm  with  the  reliability  characteristics 
of  the  Gauss-Newton  algorithm. 

4.  Develop  a  method  of  selecting  the  appropriate  design  variable  step  sizes  that  is 
more  efficient  than  the  standard  sensitivity  study. 

5 .  Develop  a  technique  that  relaxes  the  requirement  that  each  CFD  simulation  must 
be  fully  convergent. 
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6.  Reduce  computational  costs  by  using  lower  order  CFD  analysis  techniques  (MOC, 
PNS,  etc.)  in  preliminary  design  phases. 

7.  Reduce  computational  costs  by  improving  the  efficiency  of  the  CFD  analysis  code. 

8.  Reduce  the  labor  required  to  perform  each  CFD  simulation. 

Items  1  and  2  represent  research  in  areas  with  tremendous  potential.  Both  of  these 
techniques  have  the  potential  to  not  only  reduce  the  computational  effort  required  to  optimize 
a  given  design,  but  they  could  also  extend  the  allowable  design  space  without  significant 
penalty.  However,  neither  optimization  method  has  been  investigated  or  applied  extensively 
to  realistic  design  problems.  Items  3  through  S  are  ideas  to  improve  the  efficiency  of  the 
subject  direct  optimization  method.  Item  6  is  a  suggested  engineering  approach  to  reduce 
computational  costs  and  total  design  time  within  the  existing  optimization  framework.  Items 
7  and  8  address  the  CFD  analysis  portion  of  the  design  process,  which  is  where  the 
preponderance  of  computational  resources  and  manpower  is  expended.  Again,  CFD  is  not 
a  mature  science  and  we  must  continue  to  explore  more  efficient  and  more  accurate  means 
of  analysis  in  addition  to  ongoing  efforts  to  extend  the  application  of  CFD  to  more 
geometrically  and  aerodynamically  complex  configurations. 
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Figure  1.  Generic  free-jet  engine/inlet  compatibility  test  configuration. 
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Figure  2.  Design  optimization  strategy. 
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Figure  3.  Inviscid  airfoil  optimization. 
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Figure  4.  Inviscid  airfoil  surface  pressure  distribution. 
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Figure  5.  Inviscid  airfoil  derivative  sensitivity. 
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Figure  6.  Inviscid  airfoil  contour  variation. 
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Figure  10.  Invisdd  airfoil  objective  function  reduction  for  specified  surface  pressure 
distribution. 
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Figure  11.  Invisdd  airfoil  design  parameter  convergence  for  specified  surface  pressure 
distribution. 
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Figure  12.  Viscous  airfoil  optimization. 
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Figure  13.  Viscous  airfoil  target  velocity  field. 
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Figure  15.  Viscous  airfoil  objective  function  reduction. 
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Figure  17.  Separated  viscous  airfoil  contour  variation. 


Figure  18.  Separated  viscous  airfoil  initial  guess  velocity  field. 
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Figure  21.  Planar  convergent/divergent  nozzle  optimization. 


Figure  22.  Planar  convergent/divergent  nozzle  centerline  Mach  number  distribution. 
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Figure  23.  Planar  convergent/divergent  nozzle  centerline  Mach  number  comparison  to  one¬ 
dimensional  theory. 


Figure  24.  Planar  convergent/divergent  nozzle  objective  function  derivative  sensitivity. 
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Distance  Across  Reference  Plane 

Figure  25.  Planar  convergent/dlvergent  nozzle  RP  Mach  number  profiles. 
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Figure  26.  Planar  convergent/divergent  nozzle  objective  function  reduction. 
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Figure  28.  Planar  convergent/divergent  nozzle  Bezier  design  parameter,  Pj,  convergence. 
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Figure  29.  Planar  convergent/divergent  nozzle  Bezier  design  parameter,  P2,  convergence. 


Figure  30.  Planar  convergent/divergent  nozzle  total  pressure  design  parameter,  P3, 
convergence. 
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Figure  31.  Planar  convergent/divergent  nozzle  total  temperature  design  parameter,  P4, 
convergence. 
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Figure  33. 3-D  convergent/ divergent  nozzle  centerline  Mach  number  distribution  comparison. 
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Figure  34.  3-D  convergent/divergent  nozzle  derivative  sensitivity. 
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Figure  35.  3-D  convergent/divergent  nozzle  RP  Mach  number  variation  (Z  =  0). 
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Figure  38.  3-D  convergent/divergent  nozzle  design  parameter  convergence. 
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Figure  39. 3-D  convergent/ divergent  nozzle  wall  contour  convergence  for  constant  ¥  plane. 
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Figure  40. 3-D  convergent/divergent  nozzle  wall  contour  convergence  for  constant  Z  plane. 
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Figure  41.  Planar  supersonic  forebody  simulator  optimization. 
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Figure  42.  Target  design  Mach  contours  in  the  forebody  simulator/inlet  region. 
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Figure  43.  Initial  guess  Mach  contours  in  the  forebody  simulator/inlet  region. 
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Figure  45,  Planar  supersonic  forebody  simulator  RP  Mach  number  profiles. 
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Figure  46.  Planar  supersonic  forebody  simulator  objective  function  reduction. 


Figure  47.  Planar  supersonic  forebody  simulator  design  parameter  convergence. 
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Figure  48.  Comparison  of  tbe  total  pressure  component  of  the  objective  function  for  a 
shortened  fore  body  simulator. 


Figure  49.  Comparison  of  the  axial  Mach  number  component  of  the  objective  function  for 
a  shortened  forebody  simulator. 
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Figure  SO.  Comparison  of  the  vertical  Mach  number  component  of  the  objective  function 
for  a  shortened  forebody  simulator. 


.  Figure  51.  Shortened  supersonic  forebody  simulator  objective  function  reduction. 
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Figure  52.  Shortened  supersonic  forebody  simulator  RP  Mach  number  profiles. 


Figure  53.  Shortened  supersonic  forebody  simulator  variation  in  variable  geometry. 
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APPENDIX  A 

GAUSS-NEWTON  ALGORITHM 

To  formulate  the  nonlinear  least-squares  minimization  algorithm,  let  the  residuals 
ri(Pi,...,pM),  i  =  be  functions  of  the  M  design  parameters,  P.  To  minimize  n,  in 

the  least-squares  sense,  values  for  the  parameters,  P,  are  found  that  minimizes 

F(P)  =  E  {y,  -  f,(P)}2  =  E  {ri(P)}2  (A-l) 

i  =  1  l-l 


where  y,  denotes  the  N  specified  reference  plane  quantities,  and  fj  denotes  the  computed  M 
quantities  for  the  associated  design  parameters,  P. 


The  Gauss-Newton  minimization  algorithm  is  derived  by  first  substituting  a  first-order 
Taylor  series  expansion  about  P°  into  Eq.  (A-l)  for  r5(P).  This  yields 


F(P)  ~ 


N 

E 

i  -  1 


n(P°)  + 


M 

E 

j  =  i 


ftK P°) 

dPj 


APi 


(A-2) 


with  APj  defined  as  the  jth  component  of  the  vector 


AP  =  P  -  P° 


(A-3) 


Minimizing  Eq*  (A-2)  requires 


*F(P)  =  0 

apk 


(A-4) 


for  each  k.  Applying  condition,  Eq.  (A-4),  to  Eq.  (A-2)  yields  the  following  system  of 
equations: 


N 

2  E 


jri(P°)  + 


M 

E 

i  =  l 


ap]  ar‘g°>  =  0 

3Pj  J  3Pk 


(A-5) 


for  k=  1,2,3,. ...M.  Rearranging  (A-5) 


i 


N 

E 

=  i 


M 

E 

j  =  l 


*ri(P°) 

apj 


fri(P°) 

apk 


3ri(P«) 

3Pk 


(A-6) 


Defining  J  as  the  Jacobian  matrix  of  the  residual  vector,  R,  with  respect  to  the  parameter 
vector,  P,  allows  the  M  equations  defined  by  Eq.  (A-6)  to  be  written  more  concisely  in  vector 
notation  as 
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JT  JAP  =  -  JTR  (A- 7) 

where  J  is  defined  by  Eg.  (5).  This  M-by-M  system  of  equations  is  then  solved  to  determine 
the  parameter  corrections  at  each  iteration. 
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APPENDIX  B 

BROYDEN'S  JACOBIAN  UPDATE  FORMULA 

Application  of  the  Gauss-Newton  optimization  algorithm  requires  solution  of  the  system 
of  equations. 


JTJAP  =  -  JTr  (B-l) 

where  R  denotes  the  residual  vector  of  length  N,  P  denotes  the  parameter  vector  of  length 
M,  AP  denotes  the  computed  change  in  parameters,  and  J  denotes  the  Jacobian  matrix  of 
R  with  respect  to  P.  Equation  (B-l)  then  determines  a  search  vector,  AP,  and  a  linear  search 
can  be  carried  out  so  that 


Ek+I  =  Ek  +  MPk  (B-2) 

where  v  is  selected  to  reduce  the  norm  of  Rk + 1  sufficiently  below  the  norm  of  Rk.  Eliminating 
the  linear  search  is  equivalent  to  specifying  v  =  1  in  Eq.  (B-2). 

A  quasi-Newton  method  is  derived  by  assuming  at  the  kth  iteration  a  current  design  point, 
Pk,  and  an  approximation,  Bk,  to  Jk  has  been  obtained.  The  Jacobian  initial  value.  Bo,  can 
be  computed  by  a  finite  difference  approximation.  The  first-order  Taylor  series  expansion 
of  Rk4-i  gives 


Ek+l  “  Ek  +  Jk  AP  k  (B-3) 

or  in  terms  of  the  Jacobian  approximation,  Bk, 


BitAPy  =  ARk  (B-4) 

A  method  is  sought  to  approximate  Bk+j  without  reapplying  a  finite  difference 
representation.  The  method  introduced  by  Broyden  (Refs.  12  and  14),  requires  that  Bk  +  i 
satisfy  the  quasi-Newton  condition, 

Bk+1APk  =  ARk  (B-5) 

which  provides  N  equations  for  the  N  x  M  unknown  elements  of  Bk+l  in  Eq.  (B -5).  This 
is  motivated  by  examining  the  behavior  when  R  is  linear.  For  linear  R  Eq.  (B-S)  constrains 
Bk+i  to  hold  exact  derivative  information  in  the  direction  of  APk.  This  is  apparent  from 
Eq.  (B-3),  which  is  exact  for  linear  R.  The  additional  N  x  (M  - 1)  equations  are  obtained  by 
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requiring  the  change  in  R  predicted  by  Bk+i  in  all  directions  orthogonal  to  APi»  to  be  the 
same  as  would  be  predicted  by  Bk.  Symbolically,  that  provides 

Bk+ iQ_i  =  Bkqi  (B-6) 

where  the  q  i  form  a  basis  in  the  (M  - 1)  dimensional  subspace  orthogonal  to  AP*.  Thus, 

APfqs  =  0  (B-7) 

for  each  of  the  M  - 1  directions,  qi.  Define  Q  to  be  a  matrix  of  vectors  so  that 

Q  =  (APk  Q  l  0_2  Q(M-i))  (B-8) 

For  convenience  select  the  length  of  qs  as 

q/  q  i  =  APjAPk  (B-9) 

Now 

QTQ  =  (APfAPk)!  (B-10) 


which  yields 


Q-‘  = 


QT 

apJap^ 


(B-ll) 


Combining  Eq.  (B-8)  with  Eq.  (B-6)  gives 


Bk+iQ  =  BkQ  +  (Bt+jAPk  -  BkAP]j)M  (B-12) 


where 


M  =  (1  0  0  ...  0)  (B-13) 

Combining  Eq.  (B-S)  with  Eq.  (B-12)  gives 

Bk  +  iQ  =  BkQ  +  (ARk  -  BkAPk)  M  (B-14) 
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Post  multiplying  by  Q  - 1  gives 


Bk+i  =  Bk  +  (ARk  ~  BkAPOMQ  1 


Substituting  from  Eq.  (B-ll)  gives 


Bk  +  i  =  Bk  + 


(ARk  -  BkAPk)MQT 
APfAPk 


which  upon  simplification  gives 


Bk+i  =  Bk  + 


(ARk  -  BkAPk)APk 
APkAPk 


(B-15) 


(B-16) 


(B-17) 
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APPENDIX  C 

QUADRATIC  INTERPOLATION  LINEAR  SEARCH 

Following  Hartley  (Ref.  13),  let  F0  denote  the  objective  function  value  for  the  initial  value 
of  the  solution  vector,  P°,  of  a  Gauss-Newton  or  quasi-Newton  iteration.  Let  Fi  denote  the 
objective  function  value  at  the  updated  value  of  the  solution  vector,  Pl,  where  P1  is  given  by 

Pi  =  P°  +  AP  (C-l) 

with  AP  denoting  the  computed  correction  to  the  solution  vector  for  the  optimization 
iteration.  Let  v  define  a  variable  so  that  v  =  0  at  P°  and  v  =  1  at  P1.  The  objective  function 
can  now  be  approximated  in  the  interval  (0,1)  by  an  interpolation  function  with  the 
optimization  value  analytically  determined.  In  particular,  if  the  objective  function  is  evaluated 
at  v  =  m,  within  the  interval,  then  F(P°  +  vAP)  can  be  approximated  by  a  quadratic 
function  of  the  form 


f(v)  =  pv2  +  qy  +  r 


(C-2) 


If  Fm  satisfies  the  criterion 


Pm  <  Fo 


(C-13) 


and 


Fm  <  F, 


(C-4) 


then  a  minimum  of  f(v)  must  exist  in  the  interval  (0,1).  Differentiating  Eq.  (C-2)  with  respect 
to  v *  and  setting  the  result  equal  to  zero  yields 


v 


(C-5) 


with  v*  denoting  the  location  of  the  minimum  of  F(v)  in  the  interval.  The  system  of  equations 


0  0 

1  1 

,  m2  m 


F0\ 

F, 

|Fm/ 
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is  solved  to  detennine  the  coefficients  p,  q,  and  r  of  the  quadratic.  Substituting  these  coefficients 
into  Eq.  (C-5)  yields 

*  _  J_  (1  -  m2)  F0  +  m2Fi  -  Fm  — 

2  (1  -  m)  F0  +  mF,  -  Fm 

The  functional  value  f(»>*)  should  then  be  a  better  approximation  to  the  minimum  than  F|, 
which  resulted  from  the  optimization  algorithm.  Because  of  the  expense  of  function  evaluations 
within  this  research,  the  described  line  search  was  only  used  when  the  optimization  algorithm 
failed  to  produce  a  reduction  in  functional  value.  Symbolically,  if  FL  >  F0,  then  the  line 
search  was  used.  A  similar  quadratic  line  search  method  can  be  derived  for  the  interval  by 
using  Fo,  F],  and  the  f '(0).  Since  f '(0)  is  computed  to  implement  the  optimization  algorithm, 
this  approach  eliminates  the  need  for  the  evaluation  of  Fm.  This  implementation  was  not 
evaluated  herein. 
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NOMENCLATURE 
B  Jacobian  approximation 

Cj  First  algebraic  test  function  convergence  criterion 

C2  Second  algebraic  test  function  convergence  criterion 

F  Least-squares  objective  function 

gf (u)  Bernstein  basis  function  defined  by  Eq.  (10) 

J  Jacobian  matrix  of  R  with  respect  to  P  defmed  by  Eq.  (5) 

M  Total  number  of  design  parameters 

M,  Mach  number  component  along  the  x  axis 

My  Mach  number  component  along  the  y  axis 

Mz  Mach  number  component  along  the  z  axis 

N  Total  number  of  residual  components 

Nr  Number  of  reference  plane  points 

Nc  Number  of  parameter  constraints 

P  Design  parameter  vector 

Pi  Individual  design  parameter 

Pt  Total  pressure 

rj  Least-squares  residual  component 

R  Residual  component  vector 

RP  Aerodynamic  optimization  reference  plane 
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t  NACA0012  airfoil  thickness  parameter  [See  Eq.  (11)] 

Tx  Total  temperature 

u  Bemstein-Bezier  interpolation  parameter  [See  Eqs.  (8)  and  (9)] 

v  Bemstein-Bezier  interpolation  parameter  [See  Eq.  (9)] 

w  Bemstein-Bezier  interpolation  parameter  [See  Eq.  (9)] 

x  NACA0012  airfoil  axial  coordinate  [See  Eq.  (11)] 

xj  Position  vectors  of  Bezier  control  points  [See  Eq.  (8)] 

xjjk  Position  vectors  of  3-D  Bezier  control  points  [See  Eq.  (9)] 

x  (u)  Bemstein-Bezier  polynomial  [See  Eq.  (8)] 

x  (u,v,w)  3-D  Bezier  polynomial  [See  Eq.  (9)] 

y(x)  Function  used  to  define  NACA0012  airfoil  [See  Eq.  (11)] 

AP  Computed  change  in  design  parameter  vector  P 
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